home *** CD-ROM | disk | FTP | other *** search
/ The 640 MEG Shareware Studio 2 / The 640 Meg Shareware Studio CD-ROM Volume II (Data Express)(1993).ISO / prog / yamp2.zip / TESTREG.CPP < prev    next >
C/C++ Source or Header  |  1992-01-16  |  5KB  |  187 lines
  1.  
  2. /*
  3. YAMP - Yet Another Matrix Program
  4. Version: 1.2
  5. Author: Mark Von Tress, Ph.D.
  6. Date: 01/23/92
  7.  
  8. Copyright(c) Mark Von Tress 1992
  9. Portions of this code are (c) 1991 by Allen I. Holub and are used by
  10. permission
  11.  
  12. DISCLAIMER: THIS PROGRAM IS PROVIDED AS IS, WITHOUT ANY
  13. WARRANTY, EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED
  14. TO FITNESS FOR A PARTICULAR PURPOSE. THE AUTHOR DISCLAIMS
  15. ALL LIABILITY FOR DIRECT OR CONSEQUENTIAL DAMAGES RESULTING
  16. FROM USE OF THIS PROGRAM.
  17.  
  18. */
  19.  
  20. #include "virt.h"
  21.  
  22. //required global declaration for the
  23. //  matrix stack object
  24.  
  25. unsigned int _stklen = STACKLENGTH;
  26. MStack *Dispatch = new MStack;
  27.  
  28.  
  29. VMatrix &getx( int N )
  30.   // create an x matrix
  31.   {
  32.       Dispatch->Inclevel();
  33.       VMatrix x, c1, x2;
  34.  
  35.       x = Fill(N,1,0);
  36.       c1 = Fill(N,1,1.0);
  37.       for( int i=1; i<=N; i++)
  38.          x.M(i,1) = (double) (i-10);
  39.       x2 = x%x;
  40.       x = Ch( c1, Ch( x, x2) );
  41.  
  42.       // push x onto the stack
  43.       Dispatch->Push(x);
  44.       // decrement the subroutine nesting level
  45.       // and return the stack top
  46.       return Dispatch->DecReturn();
  47.   }
  48.  
  49. VMatrix &gety( VMatrix &x, VMatrix &beta)
  50.   // create a y matrix
  51.   {
  52.       Dispatch->Inclevel();
  53.       VMatrix y;
  54.  
  55.       y = x*beta;
  56.       srand(123);
  57.       for(int i=1; i<=y.r; i++) {
  58.         // use sum of 3 uniforms for an approximate
  59.         // normal random variable
  60.         int u = random(100)+random(100)+random(100)+3;
  61.         y.M(i,1) = y.m(i,1) + ((double) (u-150))/300.0;
  62.       }
  63.       Dispatch->Push(y);
  64.       return Dispatch->DecReturn();
  65.   }
  66.  
  67. VMatrix ®ression( VMatrix& x, VMatrix& y )
  68. // do a multiple linear regression
  69.       {
  70.       Dispatch->Inclevel();
  71.       VMatrix yx, reg, betahat;
  72.       int N=x.r, p=x.c;
  73.  
  74.       // solve for regression parameters using sweep
  75.       yx = Ch(y,x);
  76.       reg = Sweep( 2,p+1, Tran(yx)*yx);
  77.       // calculate mean square error
  78.       reg.M(1,1) = reg.m(1,1)/((double )( N-p));
  79.       reg.DisplayMat();
  80.  
  81.  
  82.       // solve regression using normal equations
  83.       betahat = Inv(Tran(x)*x)*Tran(x)*y;
  84.  
  85.       Dispatch->Push( betahat );
  86.       return Dispatch->DecReturn();
  87. }
  88. VMatrix &QRregression( VMatrix &x, VMatrix &y)
  89.   {
  90.      // use QR decomposition to solve regression
  91.  
  92.      Dispatch->Inclevel();
  93.      VMatrix Q, R, betahat;
  94.  
  95.      QR( x, Q, R);
  96.      betahat = Inv( Tran(R)*R )*Tran(R)*Tran(Q)*y;
  97.  
  98.      Dispatch->Push( betahat );
  99.      return Dispatch->DecReturn();
  100.   }
  101. VMatrix &GinvRegression( VMatrix &x, VMatrix &y )
  102.   {
  103.      // use Ginv to solve normal equations
  104.      Dispatch->Inclevel();
  105.  
  106.      VMatrix betahat = Ginv( Tran(x)*x )*Tran(x)*y;
  107.  
  108.      Dispatch->Push( betahat );
  109.      return Dispatch->DecReturn();
  110.   }
  111.  
  112.  
  113. VMatrix &SVDregression( VMatrix &x, VMatrix &y)
  114.   {
  115.      // use SVD to solve regression x = S Diag(V) Tran(D)
  116.      Dispatch->Inclevel();
  117.      VMatrix S, V, D, betahat, t;
  118.  
  119.  
  120.      Svd( x, S, V, D);
  121.      t = Fill( V.r, V.r, 0.0);
  122.      for( int i=1; i<=V.r; i++ ){
  123.        double tt = V.m(i,1);
  124.        t.M(i,i) = (fabs(tt) <= 1.0e-10 ) ? 0.0 : 1.0/tt;
  125.      }
  126.      betahat = D*t*Tran(S)*y;
  127.      Dispatch->Push( betahat );
  128.      return Dispatch->DecReturn();
  129.   }
  130. VMatrix &GetSerialCovar( VMatrix &R )
  131.   {
  132.     Dispatch->Inclevel();
  133.     VMatrix centered, z, spectdens, covar;
  134.     double n = (double) R.r;
  135.     centered = R - Sum( R/n ).m(1,1);
  136.     z = Fft( centered );
  137.     spectdens = Sum( z%z/n, COLUMNS);
  138.     covar = Fft( spectdens, FFALSE );
  139.     Dispatch->Push( covar );
  140.     return Dispatch->DecReturn();
  141.   }
  142.  
  143. main()
  144. {
  145.       Dispatch->Inclevel();
  146.       VMatrix x, beta("beta",3,1), y, betahat, resids, serial;
  147.       Setwid(15);
  148.       Setdec(10);
  149.  
  150.  
  151.       beta.M(1,1) = 1;
  152.       beta.M(2,1) = 0.5;
  153.       beta.M(3,1) = 0.25;
  154.  
  155.       // Also try this with 200, and note how poorly the
  156.       // ginv regression does because of colinearities
  157.  
  158.       x = getx( 55 );
  159.       y = gety(x,beta);
  160.  
  161.       betahat = regression(x,y);
  162.       betahat.Nameit( "Text book betahat");
  163.       (Tran(beta)).DisplayMat();
  164.       (Tran(betahat)).DisplayMat();
  165.  
  166.       betahat = QRregression( x,y );
  167.       betahat.Nameit( "QR betahat");
  168.       (Tran(beta)).DisplayMat();
  169.       (Tran(betahat)).DisplayMat();
  170.  
  171.       betahat = GinvRegression( x, y);
  172.       betahat.Nameit( "Ginv betahat");
  173.       (Tran(beta)).DisplayMat();
  174.       (Tran(betahat)).DisplayMat();
  175.  
  176.       betahat = SVDregression( x,y);
  177.       betahat.Nameit("SVD regression");
  178.       (Tran(beta)).DisplayMat();
  179.       (Tran(betahat)).DisplayMat();
  180.  
  181.       resids = y - x*betahat;
  182.       serial = GetSerialCovar( resids );
  183.       (Tran( Submat( serial, 5,1 ) )).DisplayMat();
  184.  
  185.    vclose();
  186. }
  187.